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O \ Abstract The basic ideas of a homotopy-based multiple-variable method is proposed 

£^' and applied to investigate the nonlinear interactions of periodic traveling waves. Math- 

r-| ' ematically, this method does not depend upon any small physical parameters at all and 

Qh. thus is more general than the traditional multiple- scale perturbation techniques. Phys- 

ically, it is found that, for a fully developed wave system, the amplitudes of all wave 
^ I components are finite even if the wave resonance condition given by Phillips (1960) 

Q\ ■ is exactly satisfied. Besides, it is revealed that there exist multiple resonant waves, 

f^ ' and that the amplitudes of resonant wave may be much smaller than those of pri- 

mary waves so that the resonant waves sometimes contain rather small part of wave 



energy. Furthermore, a wave resonance condition for arbitrary numbers of traveling 



O \ waves with large wave amplitudes is given, which logically contains Phillips' four-wave 



resonance condition but opens a way to investigate the strongly nonlinear interaction 
of more than four traveling waves with large amplitudes. This work also illustrates 
that the homotopy multiple-variable method is helpful to gain solutions with important 
physical meanings of nonlinear problems, if the multiple variables are properly defined 



5^ , with clear physical meanings. 
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1 Introduction 

In his pioneering work about nonlinear interaction of four gravity waves in deep water, 
Phillips p] gave the criterion condition for wave resonance 

ki ± k2 ± ks ± k4 = 0, (Ji ± a2 ± (73 ± (T4 = 0, (1) 

where o"j = ^/gki with ki = |kj| [i = 1,2,3,4) is the wave angular frequency, and 
kj denotes the wave number. Note that, Phillips resonance condition ([T]) works only 



for weakly nonlinear waves with small amplitudes, because ai = \fgki is the angular 
frequency of linear theory for a single gravity wave with small amplitude. 

Especially, in case of k4 = ki, Phillips [1] showed that, if ks = 2ki — k2 and 
o"3 = 2o"i — o"2, then a steady-state solution for the triad did not exist, and the ampli- 
tude of the third wave, if initially zero, would grow linearly in time. This conclusion 
was confirmed by Longuet-Higgins [2] via perturbation theory and supported by some 
experiments [3lll]. Besides, Benney [5] solved the equations governing the time de- 
pendence of the resonant modes and studied the energy-sharing mechanism involved. 

Although half century passed since Phillips' pioneering work [1] , there exist still 
some open questions about nonlinear interaction between gravity waves. First, Brether- 
ton [6] pointed out that the perturbation scheme used by Phillips breaks down for 
large time, and suggested by investigating a one- dimensional dispersive wave model 
that the amplitude of each wave component should be bounded. Are the amplitudes 
of each component waves bounded when the resonance condition is exactly satisfied? 
Besides, based on perturbation method, Phillips pQ gave the condition ([T]) of wave res- 
onance for only four waves. What is the resonance condition for more waves with large 
amplitude? It seems difficult to apply the perturbation scheme used by Phillips |T] 
and Longuet-Higgins |2] to answer these questions, because the related algebra was 
daunting and "extraordinarily tedious" , as mentioned by Phillips [7] . 

Perturbation techniques are powerful analytic tools for nonlinear equations, espe- 
cially when there indeed exists a small physical parameter, i.e. perturbation quantity: 
the solution is often expressed in a series of the small physical parameter, and the orig- 
inal nonlinear equation is transformed into a sequence of (mostly linear) sub-problems. 
Frankly speaking, perturbation techniques greatly enrich our knowledge and deepen 
our understandings about many nonlinear problems. It has a golden times before the 
times of computer when people calculated mostly by hand and wrote the "analytic" 
results on one page (or a few pages) of paper. However, some restrictions of pertur- 
bation techniques are well-known. First of all, it depends too strongly on the small 
physical parameters, but unfortunately many nonlinear problems have no such kind of 
small parameters. Besides, perturbation approximations often break down when the 
so-called perturbation quantity increases. All of these disadvantages greatly restrict 
the applications of perturbation techniques. 

In history, the calculating tools have a deep impression on the mathematical 
methods. One example is the development of numerical methods such as the finite 
element method, the finite difference method, and so on. Fortunately, we are now in 
the times of computer: a computer can do more than 10^^ fundamental operations 
and save a great lot of data in diskette within a seconds! Besides, there are some 
powerful symbolic computation software such as Maple, Mathematica, MathLab and 
so on. Today, using these software on a laptop, one can deduce lengthy formulas, 
calculate a rather complicated expression, and save analytic results on diskette in a 
seconds which might be hundreds of papers long if printed out. So, when the pencil 
and paper are replaced by the keyboard, diskette and CPU of a personal computer, 
the revolution of analytic approximation methods comes, although it is much later 



than the revolution of numerical approximation methods. 

What kind of analytic approximation methods should we have in the times of 
computer? Can we overcome the restrictions of the traditional perturbation meth- 
ods? The answer is rather positive and stimulant. Note that, the traditional concept 
of "analytic" solution came into being in the times of pen and paper, and it was tra- 
ditionally believed that a "analytic" result must be short and should be written on 
less than one paper. Today, it needs only a few seconds to save a rather complicated 
formula in diskette which might be even more than hundreds of papers. Besides, one 
needs only a few seconds to calculate such a lengthy expression by a personal com- 
puter. In fact, one spends much less time to save and calculate such lengthy formula 
by computer than to calculate and write a simple formula with half-page length by 
hand. Therefore, in the times of computer, the concept of "analytic" result must be 
modified, and an "analytic" formula can be very lengthy. Secondly, due to the high 
performance of symbolic computation of a computer, it should be easy to obtain high- 
order approximations by means of the new analytic methods. And more importantly, 
some efficient approaches should be provided to ensure the high accuracy of analytic 
approximations even for strongly nonlinear problems. In a short, the new analytic 
method should belong to the times of computer, and besides should overcome the 
restrictions of traditional methods mentioned above. 

Such a kind of analytic method has been developed in 1990s when the early 
symbolic computation software just appeared: based on the homotopy, a fundamental 
concept in topology, the homotopy analysis method (HAM) was first proposed by 
Liao [SHTG]. Different from perturbation methods, the HAM has nothing to do with 
any physical parameters. More importantly, different from other previous analytic 
methods, the HAM provides us a simple way to ensure the convergence of solution 
series. Besides, it has been proved that the HAM logically contains the previous non- 
perturbation methods such as Lyapunov artificial small parameter method, Adomian 
decomposition method, the 5-expansion method and so on. Therefore, the HAM is 
valid for strongly nonlinear problems and thus is more general. The HAM has been 
successfully applied to solve different types of nonlinear ODEs and PDEs in science 
and finance. 

The idea of multiple-scales of perturbation methods has clear physical mean- 
ings. In this article, based on the homotopy analysis method (HAM), we propose a 
multiple-variable technique for general nonlinear differential equations, which keeps 
the clear physical meaning of multiple-scales but abandons completely the small physi- 
cal parameters of perturbation techniques. In §2, the nonlinear interaction of primary 
periodic traveling waves in deep water is used as an example to describe the basic 
ideas of this approach. In §3, we give convergent series solution of a fully developed 
system of two primary waves even when Phillips' resonance condition is exactly satis- 
fied. Besides, some interesting results related to multiple resonant waves are reported. 
Especially, we reveal that Phillips resonance wave condition is mathematically equiv- 
alent to the zero eigenvalue of the eigenfunction related to the resonant wave. In 
§4.1, the expression of the eigenvalue of a fully developed wave system for arbitrary 
number of traveling waves is derived, and then a resonance condition for arbitrary 



number of small-amplitude waves is given. In §4.2, a more general resonance con- 
dition is further obtained from the physical points of view, which logically contains 
Phillips' resonance condition ([T]) but works for arbitrary number of traveling waves 
with large amplitudes. In §5, the concluding remarks, open questions and discussions 
are given. The detailed mathematical derivation is given in Appendix A, and the 
proof of a convergence theorem is given briefly in Appendix B. 



2 Mathematical description 

Let z denote the vertical co-ordinate, x, y the horizontal co-ordinates, t the time, 
z = C,{x^ y, t) the free surface, respectively. The three axises of x, y, z are perpendicular 
to each other, with the unit vector i, j, k, respectively, i.e. i-j = i-k = j-k = 0, where 
• is the multiplication dot. Assume that the vorticity is negligible and there exists a 
potential '^[x^ y, z^ t) for the velocity u that u = Vip and 

VV = 0, z<Cix,y,t) (2) 

in an incompressible flow, where 

^ .9 . d ,5 

ox oy oz 

On the free surface z = ({x, y, t), the pressure is constant, which gives from Bernoulli's 
equation the dynamic boundary condition 

^^^^^^'''"°' on'^ = C(a;,y,^), (3) 

where g is the acceleration due to gravity. Besides, z — ( vanishes following a particle, 
which gives the kinematic boundary condition 

dt dz \dx dx dy dy J 
Combining the above two equations gives the boundary condition 

+ 5^ + ^5r^ + u-VUu' =0, onz = ax,y,t), (5) 



dt^ ""dz dt V2 , 

where u = V(p and u^ = V(f ■ V(f. On the bottom, it holds 

—- = 0, as z — 7> — oo. (6) 

oz 

For details, please refer to Phillips [1] and Longuet-Higgins [2] . 

Without loss of generality, let us first consider a wave system basically composed of 
two trains of primary traveling gravity waves in deep water with wave numbers ki , k2 



and the corresponding angular frequencies 0"i,(T2, respectively, where ki x k2 7^ 
(i.e. the two traveling waves are not collinear). Due to nonlinear interaction, this 
wave system contains an infinite number of wave components with the corresponding 
wave number mki + nk2, where m, n are integers. Assume that the wave system is in 
equilibrium so that each wave amplitude is constant, i.e. independent of the time. Let 
ai,a2 denote the angles between the positive x-axis i and the wave number vectors 
ki and k2, respectively, where ki ■ k = k2 ■ k = 0, i.e. the 2;- axis is perpendicular to 
the wave numbers ki,k2. Then, 

ki = ki (cosai i + sinai j) , (7) 

k2 = k2 (cos 0:2 i + sin 0:2 j) , (8) 

where ki = |ki| and k2 = |k2|. 

Write r = xi + yj. According to linear gravity wave theory, the two trains of pri- 
mary waves traveling with the wave numbers ki and k2 are given by ai cos^i, 02 cos .^2, 
respectively, where 

^1 = ki ■ r — (Ti t = fci cos ai X + ki sin ai y — cti t (9) 

^2 = k2 ■ r — 0-2 t = k2 cos a2 x + k2 sin a2 y — cr2 t. (10) 

So, the above two variables have very clear physical meaning: for the considered 
problem, the wave profile must be a periodic function of .^1 and ^2- Mathematically, 
using these two variables, the time t should not appear explicitly for a fully develop 
wave system. In other words, one can express the potential function ip{x,y,z,t) = 
0(^5 ^2, z) and the wave surface ({x, y, t) = "^(^i, ,^2) for two trains of primary traveling 
waves. In this way, the mathematical expressions of the unknown potential function 
and wave surface are clear, with clear physical meanings. 

Obviously, it holds 

dip d(f) d^i d(f) d^2 ^■^ ^■^ 



dx d^, dx ' 96 dx ^'^ ^°^"^) W. + ^^' ^°^"^) ^' 



and similarly 



-- = (fcismai)-— + (/cssmas)^;-, 
oy ail 8^2 

dip (90 

dz dz ' 

dip d(f) d(j) 



Then, 



V(/9 



where 



Thus, 



.dtp ,d(p dip 

dx dy dz 

1 , . . .X "90 , , . . ., (90 , (90 

fci (cos ax 1 + sm a^ j ) ttt- + ^2 (cos c^2 i + sm a2 j ) tt:^ + k-- 

(94i (9^2 "9-2 

(90 (90 (90 ^ 

kiT^ + ^21^ + k-- = V0 = u, 
•941 <942 oz 



^ 1 ^ , "9 , (9 
(94i (9^2 <9z 



:ii) 



:i2) 



u^ = Vv? ■ Vv? = V0 ■ V0 



2 ( ^ 

56 



+ 



dz 



(13) 



where ki ■ k = k2 ■ k = is used, and 

ki ■ k2 = kik2 cos(q;i — 0:2). 
In general, it holds 



V0-V^ = kl 



dip ^ ^ , f d(f) dip ^ dip d(f)\ ^ J 2 ^4> 9i/j dcpdip 

^ 9.^2 ^^'^2 dz dz 



+ kl ■ k2 



+ 



for arbitrary functions 0(6, 65-2) and -0(6; 6, -2)- 
Similarly, we have 



+ kl 



(14) 



o2 J 

VV = V'0 = A;? -I + 2ki • k2 



92 



(9^0 (92 



56^6 



"2 Qt2 



56' dz 



2 ■ 



(15) 



Then, the governing equation reads 



(9^0 (9V 2 520 520 

5eF + '^^'^^ 5656+ '5^1+5;^ 



0, -oo<2;<r7(a,6), (16) 



which has the general solution 

[A cos(mei + <2) + 5 sin(ma + n6)] gl^^^+^^^l", 

where m, n are integers and A, B are integral constants. 
For the sake of simplicity, define 



(17) 



/ 



1 . . U2 

2V0-V0=- 



'^'tj +^'^-'^'5^51 + ^2 



2 ( ^ 
56 



\( 



\dz 



Using the new variables ^i and ^2, the dynamic boundary condition ([3]) becomes 

^ f d<l) d<l> \ 

On the free surface z = i]{^i,C,2), the kinematic boundary condition (jS]) reads 
2 ^V d^<P 2-9^0 d(f) f df 9/ \ ^ , ^ , 



where 
and 

On the bottom, it holds 

-^ = 0, as z ^ -00. (21) 

oz 

Given two angular frequencies cxi and (72, our aim is to find out the corresponding 
unknown potential function 0(6,6)^) ^-iid the unknown free surface 77(6,6); which 
are governed by the linear partial differential equation flT6l) subject to two nonlinear 
boundary conditions flT9l) and fl20|) on the unknown free surface z = ri{^i,C,2), and one 
linear boundary condition ( l2T]) on the bottom. Here, it should be emphasized that, by 
means of the two independent variables 6 and 6; the time t does not appear explicitly 
in the unknown potential function and wave surface. More importantly, these new 
variables have clear physical meanings. This greatly simplifies the problem solving, 
as shown later in this article. 



3 Homotopy-based approach 

As mentioned before, the two variables 6 and 6 have clear physical meanings and 
the solutions of considered problem should be periodic functions of 6 and 6- From 
physical points of view, it is clear that the wave surface should be in the form 

+00 +00 
^(6,6) = X^ X^ a,„,„cos(m6 +<2), (22) 

m=0 n=— 00 

where am,n is the amplitude of wave component cos(m6 + ^6)- Note that ( fTTj) is 
the general solution of the governing equation ( iT6l) . So, the corresponding potential 
function should be in the form 

+00 +00 



0(^1, 6, ^) = J2 Yl ^™." ^m,n(6, 6, Z) (23) 



m=0 n=— 00 



where 

^„.,„(ei, 6, z) = sin(m6 + <2) 6''"'^^+"'^^^ (24) 

and bm,n is unknown coefficient independent of ^i, ,^2; z. Note that the potential func- 
tion 0(^1,^2,^) defined by (!23|) automatically satisfies the governing equation (TT6|) . 
The above expressions are called the solution expressions of t] and 0, respectively, 
which have important role in the frame of the homotopy analysis method. 

For simplicity, define a nonlinear operator 

where the angular frequencies o"i,o"2 are given. This nonlinear operator is based on 
the nonlinear boundary condition fl20|) . Note that it contains a linear operator 



^0 (0) = -1 ^ + 2-1-2 5^^ + -2 ^ + ^^. (26) 

Let C denote an auxiliary linear differential operator with the property £[0] =0. As 
mentioned in many articles, one of advantages of the homotopy analysis method is 
the freedom on the choice of the auxiliary linear operator. Based on the results of 
linear wave theory, i.e. 



CTi ^ \/gki = d-i, 0-2 ^ \/gk2 = a2, (27) 

we choose such an auxiliary linear operator 

2 ^V 5^0 2 ^V d(f) , , 

Then, let g G [0, 1] denote the embedding parameter, cq 7^ an auxiliary param- 
eter (called convergence-control parameter), 4>oiCiyC2,z) an initial approximation of 
the potential function with V0o = 0, respectively. We construct such a parameterized 
family of equations (called the zeroth-order deformation equation) in the embedding 
parameter q E [0, 1]: 

V'0(ei,6,^;g) = o, -oo<z<r/(a,6;g), (29) 

subject to the two boundary conditions on z = 17(^1,^2; l)'- 

(1 - g) £ [0(6, 6, z; q) - Uiu 6, z)] = q Co U [0(^1, 6, z; q)] , (30) 

and 

(i-g)^(ei,6;g) 



gco<{?7(6,6;g) - - 
9 



d4>{Ci,^2,z;q) d4){Ci,C2,z;q) , 
"^1 57 ^ "^2 ^7 / 



(31) 



where 

/ = ^^k^U 6, z; q) ■ V0(6, 6, z; q). 



Besides, at the bottom, it holds 

dz 



0, as 2; — )■ —00. (32) 



Thus, using the property £ (0) = of the auxihary hnear operator fl28|) . we have 
when g = the following relationships 

0(ei,6,^;O) = 00(^1,6,^), (33) 

and 

^(6,6;0) = o, (34) 

which provide us the initial approximations of the potential function 0(^1,^2,2;) and 
the free surface ^(G)^2)- Since cq 7^ 0, when g = 1, Eqs. (l29l) to ( l32l) are equivalent 
to the original equations ( IT6l) . (iT9ll . ( I20l) and (1211) . respectively, provided 



0(6,6,^;!) = 0(6,6,^), ^(ei,6;i) = r/(6,6)- (35) 

So, as the embedding parameter q G [0, 1] increases from to 1, (j){ii^ ^2, z; q) deforms 
(or varies) continuously from the initial approximation 0o(G 5^2,-2) to the unknown 
potential function 0(^i,.^2, z), so does 77(^1,^2; q) from to the unknown wave profile 
viii^i^), respectively. Mathematically speaking, Eqs. fl29l) to fl32l) define two homo- 
topies: 

0(6, 6,^; 9) := 0o(6,6,2;) ~0(6,6,2;), 
^(6,6;?) := 0^77(^1,^2), 

This is exactly the reason why Eqs. (!29|) to (l32l) are called t/ie zeroth-order deformation 
equations. 

Note that we have freedom to choose the value of the convergence-control param- 
eter Cq. Assuming that cq is so properly chosen that the Taylor series 

+00 



te,6,^;g) = 0o(ei,6,2:) + 5^0n(ei,6,^)g", (36) 

n=l 
+00 

^(6,6;g) = E^«(^i'^2)g", (37) 

n=l 

exist and converge at g = 1, then we have due to (l35l) the homotopy- series solution 

+00 
0(6,6,2;) = </)o(6,6,^) + I]0n(6,6,^), (38) 

n=l 

+00 

^(6,6) = E^-(^i'^2), (39) 



n=l 
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where 






I9>(ei,6,2:;g) 






g=0 



ra! 



dq^' 



q=0 



are called liomotopy-derivatives. In the above formulas, the relationships ( 133|) and 
(!M|) are used. At the ?7ith-order of approximations, we have 

m 
m 



n=l 



n=l 



The equations for (f)n{C,i,^2, z) and ?7„(^i,.C2) can be derived directly from the 
zeroth-order deformation equations fl29|) to fl32|) . Substituting the series fl36l) into the 
governing equation fl29l) and the boundary condition fl32|) at bottom, and equating 
the like-power of the embedding parameter q, we have 



subject to the boundary condition at bottom 



dz 



0, as z — )■ — oo, 



(40) 



(41) 



where m > 1. It should be emphasized that the two boundary conditions ( 130|) and ( 1311) 
are satisfied on the unknown boundary ^ = 57(^1, G; <?), which itself is now dependent 
upon the embedding parameter q, too. So, it is relatively more complicated to deduce 
the corresponding equations. Briefly speaking, substituting the series (J36ll and (!37j) 
into the boundary condition (l30ll and (l3Tll with z = 57(G) ^2; Q'), then equating the 
like-power of q, we have two linear boundary conditions on z = 0: 



£ 







Co A^_^ + Xm 5*^-1 - Sm, rn> 1, 



and 



where 



and 



?/m(6,6) 



Co A^_i + Xm Vm-1, 



9 



C\ 



a 



2 ^Vm 



+ 2 0-1(72 






+ a. 



2^Vm 
9^ 



+ 9 



dz 



(42) 

(43) 
(44) 
(45) 



2=0 
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The detailed derivation of the above equations and the definitions of A^_^, S'm-i, 
Sm, Xm, Tm-i.O; 0m-i; ^m-1 ^^^ given in Appendix A. Note that, the sub-problems 
for (j)m and rjm are not only linear but also decoupled: given 4>m-i and r/m-i, it is 
straightforward to get rjm directly, and then (pm is obtained by solving the linear 
Laplace equation ( l40l) with two linear boundary conditions ( |4T1) and f j42l) . Thus, 
the high-order deformation equations can be easily solved by means of the symbolic 
computation software. 

Liao [8j proved in general that the homotopy-series solutions satisfy the original 
nonlinear equations as long as they are convergent. Similarly, we have such a theorem: 



Convergence Theorem The homotopy-series solution (|23J and ^3^ satisfy the 
original governing equation /[T6\) and the boundary conditions [T^). (d^ and /[^) . 
provided that 



-oo 



5^ At = 0, Y.^'rn = 0, (46) 

m=0 m=0 

where A^, A^ are defined by U26\) and (^^, respectively. 

A mathematical proof of the above convergence theorem is given briefly in Ap- 
pendix B. Because the terms A^ and A^ are by-products in solving high-order de- 
formation equations, the above theorem provides us a convenient way to check the 
convergence and accuracy of the homotopy-series solution. For this reason, we define 
the residual error squares 

-| pTv PIT I vri \ 

^t = -\ / E^" ^^i^^2' (47) 

£1 = -^i I U^^" ^^^^^^' (^^) 



for the mth-order approximations of (p and rj. According to this convergence theo- 
rem, the homotopy-series solution fl38|) and fl39l) satisfy the original equation and all 
boundary conditions if £^ — )■ and £^.^ — )■ as m — )■ +oo. Besides, the values of £^ 
and S^ indicate the accuracy of the mth-order approximation of and t], respectively. 



Note that the auxiliary linear operator ( l28l) has the property 

>C^m,n = [g\mki + r2k2| - (md-1 + na2f] ^^.n, (49) 



where '^rn,n is defined by ( l24l) . Thus, mathematically speaking, C has an infinite 
number of eigenf unctions \l'm.n with the corresponding eigenvalue 

Xm,n = glmki + nk2\ - {md-2 + na2) ■ (50) 

In short. 



12 



Besides, \l/m,n defined by (!24|) automatically satisfies the governing equation (TT6l) and 
the boundary condition (I2T1) at the bottom. Therefore, the inverse operator £^^ is 
defined by 

1^-' (^,n,n) = 1^, A^,„ ^ 0. (51) 

Note that the inverse operator L~^ has definition only for non-zero eigenvalue Am,n- 
When Am,n = 0, we have 

5f|mki + nk2| = (m(Ti + no'if' ^ (52) 

corresponding to the criterion of the so-called "wave resonance" mentioned by Phillips 
and Longuet-Higgins. Thus, the number of eigenf unctions with zero eigenvalue, de- 
noted by A^A, is the key for solving the problem. 

When n = 0, it holds 

Am,o = k\ {\fn\ -rr?) , 

which equals to zero only when |m| = 1 (note that m = n = corresponds to = 
and thus is not considered here). Similarly, when m = 0, the eigenvalue Ao,n equals to 
zero only when |?t,| = 1. So, there exist at least two eigenf unctions \l/i_o = e^^^siui^i 
and \l/o,i = e^2^sin^2 whose eigenvalues are zero, i.e. 

C (Ci e'''' sin ^1 + C2 e''^' sin ^2) = (53) 

for any constants Ci and C2 independent of ^1,^2 and z. Therefore, it holds N\ > 2 
in case of two primary waves. As mentioned by Phillips [1] and Longuet-Higgins p] , 
the criterion (152|) of wave resonance can be satisfied for some special wave numbers 
and angular frequencies. Thus, when ( 15^ is satisfied in case oi m = m' and n = n', 
where m' and n' are integers with m''^ + n'"^ 7^ 1, there exist three eigenf unctions 
^1,0 = e'^i'^sin^i, \l/o,i = e''2^sin^2 and 

■^m',n' = el"^''^i+"'''^l^sin(m'ei + ri'6) 
whose eigenvalues are zero, i.e. 



£ 



Ci e^^" sin ei + C2 e^2" sin ^ + C3 e''"''^^+"''^^l" sin(m'ei + n'^)! = (54) 



for any constants Ci, C2 and C3 independent of ^1, ^2 and z. Without loss of generality, 
Longuet-Higgins p] discussed a special case m' = 2 and n' = —1, corresponding to 
the eigenfunction 

^2,-1 = exp(|2ki - k2|z) sin(2ei - 6)- 

In §3 of this article, it implies m' = 2 and n' = —1 when A^^ = 3 in case of two 
primary waves, if not explicitly mentioned. 

According to the definitions ( 125]) and (H^ . it holds 



£0= (£0)1 



z=0 ■ 
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Thus, 

which gives the definition of the hnear inverse operator 



C ^ [sin(m^i + 71^2)] 



Xr 



'^m,n 7^ U. 



(55) 



Using this inverse operator, it is easy to solve the hnear Laplace equation (l40l) with two 
linear boundary conditions (14T!) and P2|) . as illustrated below. Here, we emphasize 
that the above inverse operator has definition only for non-zero eigenvalue Xm,n 7^ 0. 



3.1 In case of A^^ = 2: non-resonant waves 

In this case, there are only two eigenfunctions 

\E'i,o = exp(A;i2;) sin^i, v&o,! = exp{k2z) sin ^2 

whose eigenvalues are zero, i.e. Ai,o = Ao,i = 0. Using these two eigenfunctions and 
according to the linear wave theory, we construct the initial approximation of the 
potential function 



^0(6, 6, z) = Ao J-^ ^1,0 + Bo J-^ *o,i, 



(56) 



where An and Bn are unknown constants. 



Using formulas mentioned above or in Appendix A, the corresponding first-order 
deformation equation about the potential function 0i(^i,^2,^) reads 



V0i(ei,6,^) = o, 

subject to the boundary condition on 2; = 0: 

C (0i) = bl'^ sin(^i) + 6°'^ sin(^2) + &}'"^ sin({i + ^2) + d]'^ sin(^i - ^2) 
+ bl^' sin(2ei + 6) + dl'^ sin(2ei - 6) 
+ bl'^ sin(ei + 26) + d'/ sin(ei - 2^), 

and the boundary condition on the bottom: 



(57) 



dz 



0, 



(5^ 



(59) 



2=0 



where £ is defined by (l45l) . and U{-' ,^1 are constants. Especially, we have 



1,0 



K 



b{ 



CqAq 

cqBo 



J 

k2 



gh - ol + gh{Aohf + 2gh{Bok2f + \ sin2(ai - ^2) 

qA^k'^k^ 
gk2 -al + gk2{Bok2f + 2gk2{AokiY + ° '^ sin2(ai - ^2) 
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Since Ai^o = Ao,i = 0, according to (l55l) . it must hold 

bl'' = b'/ = 

so as to avoid the so-called "secular" terms ^i sin^i and ^2 sin ^2- This provides us the 
following algebraic equations 



{Aohf + 



2 + — ^sin^(ai -0^2) 
ZK2 



{Bok2 



fC2 . 2 / 

2 + — - sm (ai 
2ki 



tt2j 



crt 



(Tn 



-1, 



{AohY + {BohY = ^ - 1 
9k2 



(60) 
(61) 



whose solutions are 



An 



Bn 




- 1 







where 



1 9 / 

2 + — - sm (ai 
Z/C2 



"2 



k 
2 + — sin^(ai -02) 



^2-2/ N 

— sm (tti - tta) 



5f/C2 



(62) 
(63) 



-1 



-1/2 



Note that Aq and Bq have multiple values: they can be either positive or negative. 

Then, by means of the linear inverse operator C^^ defined by (l55l) . the common 
solution of 01 (^1, ^2, z) reads 



>l„/^*.o 



I J,2,l 

+ b{ 



^ 



2,1 



A 



2,1 



Bi 




^*o.H-.^' 



\&5 



A2 



&i 



1,2 



^1,1 

Ai,i 
^1,2 

Al,2 



rfi 



1,1 



di 



1,2 



^1,-1 

^1,-1 

^^1,-2 
Al -2 



(64) 



where Ai and Bi are unknown coefficients, the eigenfunction \E'm,n and eigenvalue 
Xm,n are defined by (12^ and ( 150|) . respectively. In other words, 0i is a sum ( or 
linear combination) of eigenfunctions. Note that 0i automatically satisfies the Laplace 
equation (HO!) and the bottom condition (HT!) for any constants Ai and Bi. On the 
other hand, given the initial guess 0o; it is straightforward to calculate 771(^1,^2) 
directly by means of the formula 



The above approach has general meaning. In a similar way, we can obtain 
Vm{^i,^2) and 0m(^i,^2,^), successively, in the order m = 1,2,3, and so on. Note 
that, the two unknown coefficients A^ and B^ {m > 1) can be determined exactly in 
the same way like Aq and Bq by means of avoiding the "secular" terms ^1 sin ^1 and 
^2 sin ^2- Note that only fundamental operations are needed in the above approach so 
that it is convenient to use symbolic computations to get high-order approximations. 
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Figure 1: Residual error square versus the convergence control param- 
eter Co in case of (^^ with k2/ki = 1. 



Without loss of generality, let us consider here such a special case of the two 
primary waves that 



0-1 



0-2 



71 11 

= = 1.0003, ai = 0, a2 = — , A;2 = -, 

2 6b 



(65) 



with different ratios of ki/k2- Here, the number 1.0003 is chosen so that the perturba- 
tion theory is valid with high accuracy and thus we can compare our results with those 
given by Phillips [Ij and Longuet-Higgins [2J, who suggested that the wave resonance 
(with amplitude growing in time) occurs if the criterion fl52|) is satisfied, i.e. 



fco 



0-1 



, - 0.8925, 




1.0585 



in the current case. As mentioned before, Nx = 3 when the above criterion is satisfied. 
To avoid this, let us first consider here the non-resonant waves with different wave 
number ki except k2/ki = 0.8925. 

First of all, we consider a special case ^2/^1 = 1- Without loss of generality, we 
take negative values of Aq and Bq given by fl62l) and fl63|) . respectively. To choose an 
optimal value of the convergence-control parameter cq so that the series solution of 
and 1] converge quickly, we plot the curves of the residual error squares S^ and S^ 
versus Cq, as shown in Fig. [T] When S^ and S^ contain the unknown convergence- 
control parameter cq, the related integrals are rather time-consuming. To avoid this, 
a discrete technique suggested by Liao [17] is used. It is found that the residual error 
square S^ decreases for —1.8 < Cq < and the optimal value of Cq is close to -1, as 
shown in Fig. [H Therefore, we choose Cq = —1, and the corresponding residual error 
squares of the two boundary conditions decrease rather quickly to the level 10~^^ at 
the lOth-order approximation, as listed in Table [H According to the Convergence 
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Table 1: Residual error square of two nonlinear boundary conditions in case of (j65 
with ^2/^1 = 1 



m 




^1 


1 


1.9 xlO-8 


5.1 xlO"^ 


3 


3.5 xlO-12 


1.2 xlO-9 


5 


2.0 xlO-i*^ 


4.3 xlO-i^ 


8 


2.8 xlO-22 


3.8x10-20 


10 


4.7 X 10-26 


6.4 X 10-24 



Theorem mentioned above, the corresponding homotopy-series (138|1 and (139|) are the 
solution of the problem. 

Let aio,ao,i and 02,-1 denote the amplitudes of wave components cos^i, cos ^2 
and cos(2^i — ^2)5 respectively. Obviously, Oi^o = cto,! in case of fci = ^2- As shown in 
Table [21 each wave component converges rather quickly, which agree well (see Table [3]) 
with those obtained by the homotopy-Pade method [HllTl], a kind of acceleration 
technique developed in the frame of the HAM. Besides, the analytic approximations 
of wave profile at ^2 = also converge quickly, as shown in Fig. |2l All of these indicate 
the validity of the analytic approach based on the HAM. 



Table 2: Analytic approximations of wave amplitude components in case of ([65 
k2/h = 1 



with 



Order of appr. 


C^1,0) cto,i 


^2-1 


1 
2 
3 

4 
6 

8 
10 


-0.022502 
-0.022846 
-0.022814 
-0.022816 
-0.022816 
-0.022816 
-0.022816 



0.00059739 
0.00057198 
0.00057208 
0.00057226 
0.00057226 
0.00057226 



Similarly, we can get convergent series solutions in case of ([6511 with different 
ratio of k2/ki, as shown in Table [H Note that, as the ratio of /C2//C1 decreases, the 
wave amplitude of the component cos(2,^i — ,^2) first increases monotonically to the 
maximum value at ^2/^1 = 0.8925 (the corresponding result is given by the approach 
described in §3.2), and then decreases monotonically. Note that the variation of ai^o 
and ao,i is not monotonic. It must be emphasized that |ai^o| and |ao,i| is much larger 
than |a2,i| for a// values of k2/ki, as shown in Table [H We will discuss this interesting 
phenomena later in details. 
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Table 3: Analytical approximations of wave amplitude components given by the [m, m] 
homotopy-Pade method in case of (jHSD with ^2/^1 = 1 



m 


ai, 0,(^0,1 


^2,1 


2 
3 

4 
5 


-0.022816 
-0.022816 
-0.022816 
-0.022816 


0.00057211 
0.00057226 
0.00057226 
0.00057226 



^3 = 




Figure 2: Wave profile at ,^2 = in case of fl65l) with ^2/^1 = 1- Sym- 
bols: 3rd-order approximation; Solid line: lOth-order approximation. 



Table 4: Wave amplitude components in case of fl65l) with different values of ^2/^1 



^2/^1 


01,0 


ao,i 


«2,1 


1.00 


-0.0228 


-0.0228 


0.00057 


0.95 


-0.0229 


-0.0218 


0.00084 


0.93 


-0.0230 


-0.0215 


0.00116 


0.92 


-0.0231 


-0.0214 


0.00150 


0.91 


-0.0231 


-0.0215 


0.00219 


0.905 


-0.0232 


-0.0216 


0.00290 


0.90 


-0.023 


-0.022 


0.00433 


0.8925 


-0.0205 


-0.0232 


0.00898 


0.88 


-0.0235 


-0.0179 


0.00307 


0.86 


-0.0221 


-0.0189 


0.00120 


0.85 


-0.0226 


-0.0187 


0.00087 


0.83 


-0.0229 


-0.0182 


0.00054 


0.80 


-0.0232 


-0.0172 


0.00033 


0.70 


-0.0239 


-0.0121 


0.00010 



3.2 In case of A^;^ = 3: resonant waves 



Now, let us further consider the case of (!65|) with a special ratio k2/ki = 0.8925 for 
two primary traveling waves with wave resonance. It was suggested first by Phillips [1] 
and then confirmed by Longuet-Higgins [2J that the so-called wave resonance occurs 
in this case so that the amplitude of wave component cos(2^i — ^2) grows in time, i.e. 
a2 -1 = a t , where a is a constant. 



In this case, the criterion (l52l) is satisfied, i.e. 



^|2ki 



(2ai 



^2?. 



According to fl65l) . we have 



ki = 0.703998, k2 = 0.628319, fcg = |2ki - k2| = 0.783981. 

Thus, according to (!50!) . we have an additional zero eigenvalue A2,-i = 
corresponding eigenfunction 



with the 



^ 



2,-1 



,|2ki-k2|2 



sin(2ei-e2). 



So, we have now three eigenfunctions ^1,0, ^0,1 and \E'2,-i whose eigenvalues are zero, 
i.e. Ai^o = 0, Ao,i = and A2,-i = 0. Here, it should be emphasized that the wave 
resonance condition given by Phillips [1] is mathematically equivalent to that the 
eigenvalue of nonlinear-interaction wave is zero. Unfortunately, the new zero eigen- 
value A2,-i breaks down the approach mentioned in §3.1, because the term ^E'2.-i/A2,-i 
in ( 1M1) becomes infinite: this is exactly the reason why Phillips [1] and Lenguet- 
Higgins [2] suggested the existence of the so-called wave-resonance with amplitude 
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growing in time, which, however, is physically impossible from the view-points of 
wave energy. 

Can we avoid such kind of wave-resonance with amplitude growing in time? 

Note that, in case of non-resonant waves investigated in §3.1, the initial guess 
( l56ll is a linear combination of the two eigenfunctions \E'i.o and \l'o,i whose eigenvalues 
Ai^o and Ao,i are zero. In the current case of wave resonance, the only difference is 
that we have an additional eigenfunction ^2,-1 whose eigenvalue A2,-i is zero, too. 
As mentioned in many other publications P4T6]. the HAM provides us with great 
freedom to choose the initial guess. With such kind of freedom, why not use these 
three eigenfunctions (with zero eigenvalue) to express the initial guess 0o? In other 
words, we can express the initial guess 0o by all eigenfunctions whose eigenvalues are 
zero, i.e. 



^0(6,6,^) = A 



^1,0 + Bo 



^0,1 + Co 



^ 



2,-1; 



(66) 



where Aq, Bq, Cq are unknown constants independent of ,^1, ^2 and z. Similarly, substi- 
tuting the above expression into the deformation equations PUI) to fH2|) . we have the 
same first-order deformation equation fl57j) with the same boundary condition fl59l) at 
bottom for 0i, but a more complicated boundary condition on ^ = 0, i.e. 



= 6}'° sin(ei) + 6?'' sin(e2) + 6f sin(2ei) 
+ 6}'^ sin(^i + ^2) + dl'^ sin(^i - ^2) 
+ bl'^ sin(2ei + 6) + ^'' sin(26 - 6) 



+ 
+ 
+ 
+ 
+ 



&rsin(36) 



1,2 



b[''sm{^i + 2^2) + d\''sm{^i~2^2 



df sin(2ei - 2^2) + df' sin(2ei - 36) 
d{ sin(3,^i 
d{ sin(46 
d{ sin(56 



6) + rfr sin(3ei - 26) 

6) + df sin(4ei - 26) + d"/ sin(46 



36 



26) + c/rsin(66-36 



(67) 



where hi' ,di' are constant coefficients, and the linear operator C is defined by fl45|) . 
Note that there exist now three zero eigenvalues, i.e. Ai.o = 0, Ao,i = and A2.-1 = 0. 
Therefore, according to the definition fl55l) of the inverse operator C~^, not only the 
two coefficients b^' , b^' but also the additional coefficient d{ must be zero. Enforcing 



ri,o 



b['^ = 0, b\ 



0,1 



0, d^{^ 



0, 



we obtain a set of nonlinear algebraic equations 

12.7576^2 + 20.3675^2 + 2,imQSCl + 25.67185oC'o _ 
24.1456^2 + 9.60045| + 30.0398C| + 14.9558A|Co/5o 
26.9621Ag + 21.61165^ + 16.6956^2 + 10.7158AgSo/Co 



0.0154469, 
0.014593, (6^ 
0.0163008 



for the special case mentioned above. The set of these nonlinear algebraic equations 
has four complex and twelve real solutions. Because the complex solutions have no 
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physical meanings, we list only its twelve real roots in Table [51 It is found that the 
twelve roots fall into three groups, and different groups give different solutions, as 
shown later. After solving this set of nonlinear algebraic equations, the initial guess 
00 is known and therefore it is straightforward to get rji directly by means of ( H3ll . 
More importantly, on the right-hand side of Eq. (1671) . the terms sin,^i, sin ,^2 and 
especially sin(2^i — ^2) disappear now. Then, using the inverse operator fl55|) . it is 
straightforward to get the common solution of the first-order approximation 



i;2,o 



^2,0 \ , 73,0 / ^3,o\ , Tl,l / ^l,l\ , Tl,l / ^1 



+ K T^ + ^1 T^ + ^1 T^ + ^i 



-^2,0 / V -^3,0 / V -^1,1 / V -^1 



-1 






-^2,-3/ V -^3,-1 / V -^3,-2 / \ A4,_i 

^ 74,2 f ^4,-2 \ ^ 74,3 f ^4,-3 \ ^ 1^,2 f ^5,-2 \ ^ 1G,3 f ^6,-3 \ ,^q^ 

+ '^ Ia^J+'^ Ia^J+'^ Ia^J+'^ Ia^J- ^''^ 

It should be emphasized that all eigenvalues Xm,n listed in the above expression are 
nonzero so that 0i is finite. More importantly, all coefficients in the above expression 
are independent of the time so that the corresponding wave profile does not grow in 
time! Note that, like the initial guess 0o defined by fl66|) . the common solution 0i 
given by ( l69|) has three unknown coefficients Ai, Bi and C*!, which can be determined 
similarly by avoiding the "secular" terms in 02- So, the above approach has general 
meanings. Therefore, in a similar way, we can get rjm and 0m successively, where 
m = 1, 2, 3 and so on. 

Note that the wave amplitude components ai o and ao,i in Table|l]are negative. To 
calculate the corresponding wave amplitude components in case of fl65l) with ^2/^1 = 
0.8925, we choose the 2nd root in Group-I, i.e. 

Ao = -0.0156112, Bo = -0.0282054, Co = 0.00849726. 

Similarly, we can choose an optimal value of the so-called convergence-control param- 
eter Co by plotting the curves of the residual error square £^ versus co, as shown in 
Fig. [3], which indicates that the series solution converges in the region —1.6 < co < 
and that the optimal value of cq is close to -1. For simplicity, we take cq = —1. The 
residual error squares £^ and £^ of the two boundary conditions decrease rapidly to 
the level 10^^^ (at the 20th-order approximation), as shown in Table [61 According 
to the Convergence Theorem proved in Appendix B, the homotopy-series (1381) and 
(J2S]) satisfy the original governing equation (1161) and all boundary conditions (I19p . 
f[20|) and f[2T|) . Besides, it is found that the corresponding wave amplitude compo- 
nents aifi, ao,i and 02,-1 converge to -0.02051, -0.023212, 0.0089752, respectively, as 
shown in Table [7l To confirm the convergence, we further employ the homotopy-Pade 
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Table 5: Roots of Eq. ([68]) 



Series number 
of roots (K) 


Ao 


Bo 


Co 


1 (Group-I) 

2 

3 

4 


-0.0156112 

-0.0156112 

0.0156112 

0.0156112 


0.0282054 
-0.0282054 

0.0282054 
-0.0282054 


-0.0084973 
0.0084973 

-0.0084973 
0.0084973 


5 (Group-II) 
6 

7 
8 


-0.0155774 

-0.0155774 

0.0155774 

0.0155774 


-0.0141927 
0.0141927 

-0.0141927 
0.0141927 


-0.0113800 
0.0113800 

-0.0113800 
0.0113800 


9 (Group-Ill) 

10 

11 

12 


-0.0155626 

-0.0155626 

0.0155626 

0.0155626 


0.0106109 
-0.0106109 

0.0106109 
-0.0106109 


-0.0226353 
0.0226353 

-0.0226353 
0.0226353 



technique [8l[Tl] to accelerate the convergence and obtain the same convergent wave 
amplitude components 

ai = -0.0205119, ao,i = -0.0232118, a2,_i = 0.0089752, 

as shown in Table [81 Furthermore, the corresponding wave profile converges quickly, 
too, as shown in Figs. [H and [51 Therefore, we indeed get convergent series solution 
of resonant waves with constant amplitudes even when the resonant condition ( 152|) is 
exactly satisfied. 

Combining the above result with those listed in Table [H we obtain the whole pat- 
tern of the dimensionless wave amplitude component A;3(a2,_i) versus k2/ki in case 
of fl65|) . as shown in Fig. [61 It is true that, as the ratio k2/ki goes to 0.8925, corre- 
sponding to the criterion fl52l) of wave resonance, the dimensionless wave amplitude 
component k^ 0-2-1 arrives its maximum. Besides, the resonant wave profile in case 
of k2/ki = 0.8925 becomes more complicated, if compared with the non-resonant one 
in case of k2/ki = 1. However, it should be emphasized that the amplitude |a2,-i| of 
the wave component cos(2^i — ^2) is a finite constant, even if the resonance condition 
(!52|) is satisfied exactly. Besides, it is surprising that, in case of k2/ki = 0.8925, the 
amplitude 1 02,-1 1 of the resonant wave component is even much sma//er than the wave 
amplitudes |ai.o| and |ao,i| of the two primary waves! 

The above results are obtained by using the 2nd root of Group I in Table [3 Sim- 
ilarly, using different roots in Table [5l we can search for the corresponding convergent 
series solutions. It is found that the four different roots of each group in Table [51 give 
the different wave amplitude components Oi^, ao,i and 02-1- But, they have the same 
absolute values |ai,o|, |oo,i| and |a2,-i|, as listed in Table [HI Considering the fact that 
the wave energy spectrum is determined by the amplitude square of wave components. 
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we regard the four different solutions in each group as the same. Thus, in case of (!65l) 
with k2/ki = 0.8925, there are three different resonant-wave patterns with different 
wave energy spectrums. The resonant wave profiles of Group II and III are as shown 
in Figs. [7] and [HI Here, we would like to emphasize that the amplitude |a2,-i| of the 
resonant wave is the smallest in Group I, and is the middle in Group II, although it is 
the largest in Group III. So, the amplitude |a2,-i| of the resonant wave is not special 
at all: it is just normal as the wave amplitude components |ai,o| and |ao,i| of the two 
primary waves. What we would like to emphasize here is that, for a fully developed 
wave system, there exist multiple solutions when the resonance condition is exactly 
satisfied. Besides, the resonant wave amplitude may be much smaller than primary 
wave amplitudes. These interesting results have not been reported, to the best of our 
knowledge. 

Let n denote the sum of amplitude square of all wave components and write 

Ho = fli^o + '^0,1 + '^2,-1- 

It is found that Uo/U = 98.82%, 98.27% and 99.76% for Group I, II and III in case 
of fl65l) with k2/ki = 0.8925 for resonant waves, respectively. This is mainly because 
amplitudes of other wave components such as ai,_2 are much smaller, as shown in 
Table O Thus, these three wave components nearly contain the whole wave energy. 
Note that, given two primary traveling waves with wave numbers ki and k2, there 
exist an infinite number of different wave components am,nCos{mC,2 + ^^^2) with the 
wave number mki + n\<.2, where m and n are arbitrary integers. Let 

K = {mki + n]<.2\m, n are integers} 

denote a set of all these wave numbers. Each wave number in K corresponds to an 
eigenfunction defined by 1^^ with an eigenvalue defined by (15U]) . Our computations 
suggest that, for a fully developed wave system with small amplitudes, the main of 
wave energy focuses on the wave components whose eigenvalues are zero (or close to 
zero). This provides us an alternative explanation for the so-called wave resonance. 
According to this explanation, the resonant wave a2._i cos(2^i — ^2) is as important 
as the two primary waves ai^ cos^i and ao,i cos^2- 

The amplitudes of the resonant waves (related to Group I) in case of cti = 0, 02 = 
it/36, k2/ki = 0.8925 with different ratios of Oxj \/gk\ = a2/\/gk2 are as shown in 
Table [TOl The corresponding wave energy distributions are given in Table [131 It 
is found that the two primary waves contain most of the wave energy in this spe- 
cial case. Besides, as the ratio of ai/y/gki = 02I \fgk2 increases, the resonant wave 
02,-1 cos(2^i — ^2) contains less and less percentage of the whole wave energy. Espe- 
cially, when Gij \fgk\ = 02! \fgk2 = 1.0008, the resonant wave of Group I contains 
only 2.22% of the whole wave energy. This result is interesting, but a little surprising, 
because the resonance wave is traditionally supposed to have a large wave amplitude 
and thus to contain the main of wave energy. We will attempt to explain this phe- 
nomena in §4. The amplitudes of the resonant waves related to Group II are listed 
in Table [11] and the corresponding wave energy distribution is given in Table [TH It 



23 




Figure 3: Residual error square versus the convergence control param- 
eter Co in case of f l65|) with ^2/^1 = 0.8925. 



is found that the resonant waves of Group II have the comparable wave amplitudes 
with the comparable percentage of wave energy to one of primary waves. Among 
three groups, there exists only one group (i.e. Group III) such that the resonant 
waves have the largest wave amplitude and besides contain the main part of wave 
energy, as shown in Table [12] and Fig. |9l Note that, the primary and resonant waves 
contain the most part of wave energy, especially when all wave amplitude components 
are very small, as shown in Fig. [101 However, as the wave amplitudes increase, the 
primary and resonant waves contain less and less percentage of wave energy, as shown 
in Fig. [TOl This indicates that Phillips' wave resonance condition ([1]) might hold only 
for small-amplitude traveling waves. 

The above results have general meaning, although they are obtained in a special 
case OnS])- These results strongly suggest that, for a fully developed system of two 
traveling waves, all wave amplitudes do not grow linearly in time even if the wave 
resonance condition is exactly satisfied. This conclusion is also true for arbitrary 
number of traveling waves, as shown below. Currently, by means of DNS (direct 
numerical simulation) of the evolution of nonlinear random water waves fields with 
a continuous spectrum, Annenkov et al [H] investigated the role of exactly resonant, 
nearly resonant and non-resonant wave interactions, and their results indicate that 
the amplitudes of wave packets tend to constants. Their results, although obtained 
for a continuous wave spectrum, support our conclusions mentioned above. 
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Table 6: Residual error square of the two boundary conditions in case of ([65 
fca/fci = 0.8925 



with 



m 






1 


4.0 xlO-^ 


4.8x10-^ 


3 


1.8x10-8 


1.0x10-6 


5 


2.4 xlO-i° 


1.2x10-8 


8 


1.2 X 10-13 


1.5x10-11 


10 


1.3 xlO-i^ 


1.2x10-12 


15 


1.2 xlO-i^ 


3.4x10-16 


20 


1.2 X 10-20 


1.9x10-1*^ 



Table 7: Analytic approximations of wave amplitude components in case of ([65 
k^/ki = 0.8925 



with 



Order of appr. 


ai,o 


ao,i 


02,-1 


1 


-0.015616 


-0.028214 


0.0084999 


3 


-0.020587 


-0.023469 


0.0094143 


5 


-0.020533 


-0.023231 


0.0090083 


7 


-0.020504 


-0.023220 


0.0089755 


9 


-0.020511 


-0.023213 


0.0089757 


11 


-0.020512 


-0.023212 


0.0089753 


13 


-0.020512 


-0.023212 


0.0089752 


15 


-0.020512 


-0.023212 


0.0089752 


18 


-0.020512 


-0.023212 


0.0089752 
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Table 8: Analytical approximations of wave amplitude components given by the [m, m] 
homotopy-Pade method in case of fl65l) with ^2/^1 = 0.8925 



m 


«1,0 


ao,i 


02,-1 


2 
3 

4 
5 
6 
7 
8 
9 


-0.0206010 
-0.0204911 
-0.0205102 
-0.0205122 
-0.0205119 
-0.0205119 
-0.0205119 
-0.0205119 


-0.0232058 
-0.0232251 
-0.0232172 
-0.0232118 
-0.0232118 
-0.0232118 
-0.0232118 
-0.0232118 


0.0089660 
0.0089568 
0.0089780 
0.0089752 
0.0089752 
0.0089752 
0.0089752 
0.0089752 



Table 9: Multiple amplitudes of resonant waves in case of f l65|) with k2/ki 
given by different roots of Eq. f l68|) listed in Table |5l 



0.8925 





ai,o 


ao,i 


02,-1 


ai-2 


Group I 
Group II 
Group III 


0.02051186921 
0.01475607438 
0.00971236473 


0.02321179687 
0.01002488089 
0.01032248128 


0.00897520547 
0.01464333477 
0.02576462018 


0.00022907754 
0.00096077706 
0.00059968416 



Table 10: Amphtudes of resonant waves ( Group I ) in case of ai 
7r/36, k2/ki = 0.8925 with different ratios of a\l ^gk\ = oij \fgk2. 



0,02 



<yiNgki,(y2/^gk2 


ki,o 


ao,i 


^2,-1 


1.0001 


0.0117 


0.0141 


0.0058 


1.0002 


0.0166 


0.0195 


0.0078 


1.0003 


0.0205 


0.0232 


0.0090 


1.0004 


0.0238 


0.0259 


0.0095 


1.0005 


0.0266 


0.0278 


0.0096 


1.0006 


0.0290 


0.0291 


0.0092 


1.0007 


0.0312 


0.0301 


0.0085 


1.0008 


0.0331 


0.0308 


0.0076 
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Table 11: Amplitudes of resonant waves ( Group II ) in case of ai 
tt/36, k2/ki = 0.8925 with different ratios of a\j \fgk\ = 02! \fgk2- 



0,02 



<^iNgki,a2/y/gk2 


ki.o 


^0,1 


^2-1 


1.0001 


0.0090 


0.0063 


0.0084 


1.0002 


0.0124 


0.0085 


0.0119 


1.0003 


0.0148 


0.0100 


0.0146 


1.0004 


0.0166 


0.0111 


0.0169 


1.0005 


0.0181 


0.0120 


0.0188 


1.0006 


0.0194 


0.0126 


0.0206 


1.0007 


0.0204 


0.0130 


0.0222 


1.0008 


0.0214 


0.0134 


0.0237 



Table 12: Amplitudes of resonant waves ( Group III ) in case of ai 
7r/36, k2/ki = 0.8925 with different ratios of Oxj ^[gkl = 02! \fgk'2- 



0,a2 



(^iNgki,a2/^gk2 


|ai,o 


ao,i 


02,-1 


1.0001 


0.0060 


0.0060 


0.0148 


1.0002 


0.0082 


0.0085 


0.0210 


1.0003 


0.0097 


0.0103 


0.0258 


1.0004 


0.0108 


0.0118 


0.0298 


1.0005 


0.0116 


0.0131 


0.0333 


1.0006 


0.0122 


0.0143 


0.0365 


1.0007 


0.0125 


0.0153 


0.0394 


1.0008 


0.0127 


0.0163 


0.0420 



Table 13: Wave energy distribution of resonant waves 
0, 0:2 = 7r/36, k2/ki = 0.8925 with different ratios of ai/ 




in case of «i 



(^i/Vgki,(T2/Vgk2 


«?.o/n 


«o,i/n 


«2.-i/n 


Ho/n 


1.0001 


36.97% 


53.95% 


8.99% 


99.91% 


1.0002 


38.40% 


52.74% 


8.42% 


99.56% 


1.0003 


39.97% 


51.19% 


7.65% 


98.82% 


1.0004 


41.59% 


49.28% 


6.70% 


97.57% 


1.0005 


43.10% 


47.09% 


5.60% 


95.78% 


1.0006 


44.38% 


44.75% 


4.43% 


93.56% 


1.0007 


45.65% 


42.73% 


3.58% 


91.95% 


1.0008 


46.23% 


39.07% 


2.22% 


87.52% 



27 



Table 14: Wave energy distribution of resonant waves ( Group II 
0, ai = 7r/36, k2/ki = 0.8925 with different ratios of oxj \fgk[ = 0-2/ 



in case of ai 




o-i/\/9h,o-2/V9h 


«lo/n 


<i/n 


al_jU 


Ho/n 


1.0001 


42.30% 


20.43% 


37.07% 


99.79% 


1.0002 


41.26% 


19.50% 


38.44% 


99.20% 


1.0003 


40.17% 


18.54% 


39.56% 


98.27% 


1.0004 


39.05% 


17.56% 


40.44% 


97.06% 


1.0005 


37.93% 


16.56% 


41.12% 


95.61% 


1.0006 


36.79% 


15.54% 


41.64% 


93.97% 


1.0007 


35.65% 


14.51% 


42.02% 


92.18% 


1.0008 


34.49% 


13.46% 


42.29% 


90.24% 



Table 15: Wave energy distribution of resonant waves ( Group III 
0, 0:2 = 7r/36, k2/ki = 0.8925 with different ratios of o\j \fgk\ = 02I 



in case of ai 




(Til\/9ki,(r2l^/gk2 


«?.o/n 


<i/n 


«2.-i/n 


Ho/n 


1.0001 


12.13% 


12.47% 


75.37% 


99.97% 


1.0002 


11.53% 


12.38% 


75.99% 


99 


90% 


1.0003 


10.88% 


12.29% 


76.58% 


99 


76% 


1.0004 


10.20% 


12.20% 


77.15% 


99 


55% 


1.0005 


9.47% 


12.11% 


77.68% 


99 


26% 


1.0006 


8.72% 


12.01% 


78.15% 


98 


88% 


1.0007 


7.94% 


11.90% 


78.53% 


98 


38% 


1.0008 


7.15% 


11.78% 


78.82% 


97 


75% 
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Figure 4: Resonant wave profile at ^2 = in case of (!65|) witli 
^2/^1 = 0.8925. Symbols: 5th-order approximation; Solid line: lOtli- 
order approximation. 
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Figure 5: Resonant wave profile at ^1 = in case of (!65|) with 
^2/^1 = 0.8925. Symbols: 5th-order approximation; Solid line: lOth- 
order approximation. 
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k/k, 



Figure 6: The dimensionless wave amplitude component ^3(02,-1) ver- 
sus k2/ki in case of fl65l) . Filled circle: result when ^2/^1 = 0.8925; 
Open circles: results when /C2//C1 7^ 0.8925. 




Figure 7: Resonant wave profile at ^2 = in case of f l65|) with ^2/^1 = 
0.8925. Solid line: the lOth-order of approximation of Group II; Dashed 
line: the lOth-order of approximation of Group III; Symbols: 6th-order 
approximations . 
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Figure 8: Resonant wave profile at ^i = in case of (16 5 p witfi A;2/^i = 
0.8925. Solid line: the lOth-order of approximation of Group II; Dashed 
line: the lOth-order of approximation of Group III; Symbols: 6th-order 
approximations . 



Group III 



Group II 




<S,/ (gk,) 



Figure 9: a^ _i/n versus Oxj \fgk[ in case of fl65|) when ^2/^1 = 0.8925 
and cf\l \fgk[ = a-^/ \Jgk2 for different group of solutions. Solid line: 
Group I; Dashed line: Group II; Dash-dotted line: Group III. 
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Group III 




Group II 



Group I 



<^,/ {gk,}" 



Figure 10: Ilo/n versus ai/y/gk[ in case of ( 165|) when /c2/^i = 0.8925 
and (7i/ yjgki = o^j \fgk2 for different group of solutions. Solid line: 
Group I; Dashed line: Group II; Dash-dotted line: Group III. 
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4 Resonance condition of arbitrary number of waves 

4.1 Resonance condition for small-amplitude waves 

In §3, we shaw that, for a fully developed system of two primary traveling waves, 
Phillips' resonance condition for small-amplitude waves is exactly equivalent to the 
zero eigenvalue of the eigenfunction related to the resonant wave. This conclusion 
has general meanings, and can be easily expanded to give a resonance condition for 
arbitrary number of travel waves with small amplitude. The key is to give an explicit 
expression of the eigenvalue in case of arbitrary number of traveling waves. 

Let us consider the nonlinear interaction of k, periodic traveling waves with small 
amplitudes, where 2 < k < +oo is an arbitrary integer. Let k„ and cr„ (1 < n < k) 
denote the given wave number and angular frequency of the nth periodic traveling 
waves in deep water. Define the variable 

which has clear physical meanings. Then, 

ip{x, y, z, t) = 0(6, 6, ■ ■ ■ , U, z), C{x, y, t) = 77(6, 6, ■ ■ ■ , Q- 
Similarly, we have 

v.- (Ek..f)-k|--v.. (n, 

\m=l / 



K hi 



\m=l n=l / ^ ' 

V7^ V7 / /^V^V^i 1 -90 (9^ \ (90(9^ 



Note that 



^ mi,m2,--- ,mfi,z CXp 

n=l 

satisfies the Laplace equation V0 = 0, i.e. 



/ ^ TTifjKjj 



z I sm 

n=l 



^ rrinin (75) 



V^^ = 



The two nonlinear boundary conditions on the free surface can be written by means of 
the operators defined above in a similar way. Then, we can construct the zeroth-order 



33 



deformation equations and the corresponding high-order deformation equations in a 
similar way as mentioned in §3. Although it seems that the governing equations and 
boundary conditions become much more complicated in form than the original ones 
by means of these variables, these multiple variables have very clear physical meanings 
which in fact greatly simplify solving the problem, as described below. 

Similarly, we choose such an auxiliary linear operator 



£0 



EE^- 



92 



92 



0"^ 



\m=l n=l 



dimdir. 



dz^ 



(76) 



where 



a„ 



a/^ 



h 

rri) '^m 



is based on the linear theory for small-amplitude waves. The above auxiliary linear 
operator satisfies 



£(^ 



mi,m2,--- ,171^ 



'^mi,m2,--- ,171,^ ^ mi,m2,--- ,mi^y 



where 



'^mi, 



m2,--- ,mK 



/ ^ ITT'Ti'^n 



n=l 






''^n^n 



(77) 



(7^ 



is the eigenvalue and \I/mi,m2,--,m« defined by ( !75l) is the eigenfunction of the linear 
operator C defined by (1761) . Similarly, the inverse operator of ( 1761) satisfies 



C-' (*„„ 



m 



1711,1712, ■■■ ,mre 



m2,--- ,mK 



A 



A 



I7ii,m2,--- ,mK 



.^0. 



(79) 



mi,m2,--- ,m,K 



A 



Note that the inverse operator (179|) has definition only for non-zero eigenvalue 
Besides, the eigenvalues of all primary traveling waves are zero, i.e. 



7n\,m2,--- ,171^' 



there exist at least n zero eigenvalues for n primary waves. Thus, the so-called 
wave resonance occurs when there are more than k zero eigenvalues. So, enforcing 



A 



mi,m2. 



gives the resonance condition 



y^ rrinkn 



n=l 



E 



,n=l 



m„,(Tr, 



E"^« > ^' 



^0) 



n=l 



where an = y/g k^ with k^ = |k„| is based on the linear theory for small-amplitude 
waves. Note that (152|) is a special case of the above resonance condition. Besides, the 
above formula contains the resonance condition given by Phillips [1] and thus is more 
general. 

Assume that, for given k, primary traveling waves, there are A^a ^ 1^ eigenfunctions 
whose eigenvalues are zero. When A^a = 1^, there is no wave resonance. However, 
when A^A > 1^, wave resonance occurs: the wave energy transfers greatly between the 
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resonant wave and primary ones. For simplicity, let \E'5^ (1 < m, < Nx) denote the 
mth eigenfunction with zero eigenvalue. According to ( 177|) . it holds 



\m=l J 



Nx>K 



for any constant Am- So, we can always choose such an initial guess that 



f82) 



where -Bo,m is unknown. Similarly, the Nx unknown constants -Bo,m (1 < ?7i < Nx) 
are determined by avoiding the "secular" terms in 0i. Besides, owing to fIST]) . the 
common solution of 0i contains Nx unknown constants Bi^ {1 < m < Nx), which 
are similarly determined by avoiding the "secular" terms in 02- In this way, one can 
solve the related high-order deformation equations successively, and an optimal value 
of convergence-control parameter cq can be chosen so as to ensure the homotopy-series 
convergent quickly. In theory, the above approach is general, and works for arbitrary 
number of primary periodic traveling waves with small amplitudes. It provides us a 
new way to investigate the weakly nonlinear interactions of more than four primary 
traveling waves with small amplitudes. 



4.2 Resonance condition for large wave-amplitude 



Note that the general wave resonance condition flHUl) holds only in case of o"„ = 
\fgKi with kn = |k„|, corresponding to small- amplitude gravity waves. What is 
the resonance condition for arbitrary number of traveling gravity waves with large- 
amplitude? 




Figure 11: The resonance of a simple pendulum 



To answer this question, we should consider the physical meanings of f l8UI) . In 
general, the so-called resonance of a dynamic system occurs when the frequency of 
an external force (or disturbance) equals to the "natural" frequency of the dynamic 
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system. For example, let us consider the resonance of a simple pendulum, as shown 
in Fig. [m where F = A cos{ujt + a) is the external force with the frequency u 
and the phase difference a. When the maximum angle of oscillation Omax is so small 
that sin 6' ~ 9, the simple pendulum has a natural frequency cuq ~ \/9/l- So, if the 
frequency u of the external force F is equal to the natural frequency uq of the simple 
pendulum, i.e. cu = \fgjl, the total energy of the pendulum (and therefore Qrnax) 
quickly increases in case of the phase difference a = (or decreases in case of a = tt): 
the so-called resonance occurs. However, cuq ~ ydl^ is only valid for small dmax'- the 
natural frequency cjq increases as Qmax becomes larger. So, as the maximum angle 
of oscillation Qmax becomes so large that the natural frequency cuq departs more and 
more from the frequency u = ygjl of the external force F, then the simple pendulum 
gains less and less energy from the external force: the maximum angle of oscillation 
Qmax stops increasing when the simple pendulum can not gain energy from F any 
more in a period of oscillation. 

The phenomenon of gravity wave resonance is physically similar to it in essence. 
For a single traveling wave with the wave number k' and the "natural" angular fre- 
quency (Jq, the resonance occurs when there exists an "external" periodic disturbance 
with the same angular frequency a', i.e. a' = o'^. It should be emphasized that this 
resonance mechanism is physically reasonable even for large wave amplitude. 

Let us consider k primary traveling waves with wave number k^ and angular 
frequency cr„, where \ <n < k. Due to nonlinear interaction, there exist a system of 
an infinite number of wave components 



cos 



^Y^mnin 



,n=l 



where m„ is an integer that can be negative, zero, or positive. Note that 



^mnin= 5Z 



TTlnk-n 



n=l 



.n=l 



So, 



k' = E 



■ r 



TTln '^n 



vn=l 



nin cr. 



n 



t. 



n=l 



is the wavenumber and 



a 



E 

n=l 



m„, ar, 



^3) 



(84) 



is the corresponding angular frequency of the nonlinear-interaction wave. For the sake 
of simplicity, we call k' the nonlinear-interaction wavenumber and a' the nonlinear- 

K 

interaction angular frequency, respectively, where ^ m'^ > 1. 

n=l 



For small- amplitude waves, we have an ~ \/g kn = cTn, which leads to 



E 

n=l 



m„. (Tr. 



E 

n=l 



m„ o.„ 



a 



^5) 
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Substituting the above expression and (183!) into (IHOj) . we have the resonance condition 
(for K small-amphtude primary waves) in the form: 



^Ik'l 



,'2 



a 



I.e. 

a' = ^WV (86) 

The above resonance condition clearly reveals the physical relationship between the 



nonlinear-interaction wavenumber k' 



^ rrtn k„ and the nonlinear-interaction an- 

n=l 



gular frequency a' 



E ^n 0-n 



n=l 



Mathematically, the wave resonance condition (1861) can be derived in the frame 
of the HAM. Note that the HAM provides us great freedom to chose the auxiliary 
linear operator. The auxiliary linear operator C defined by fl76|) contains the term 
(Ti = y^g\k~\ which has physical meaning only for small-amplitude waves. So, for 
large-amplitude waves, we should replace the term a^ = ^/^jk^ in fl76|) by the given 
angular frequency ctj. In other words, for primary waves with large-amplitudes, we 
should choose the auxiliary linear operator 



£0 



EE 

<m=l n=l 



d^ 



CTm. CTr. 



d^mdCn 



9 



92 



dz^ 



(87) 



where cxj is the angular frequency of the ith primary wave. In fact, the above linear 
operator comes from the linear part of the nonlinear boundary condition ([5]). The 
eigenvalue of the above linear operator reads 



A 



mi,m2,--- ,"Ik 



/ ^ l^n'^n 



n=l 



E 



,n=l 



m„cr, 



n'^n 



Enforcing the above formula to be zero, we obtain 



E 

n=l 



ninan 



\ 



/ J ^T^n'^n 



n=l 



^9) 



which is exactly the wave resonance condition a' = \^gk.' defined by ( l86l) . 

Let ctq denote the "natural" angular frequency of a single traveling wave with 
the wavenumber k' and the wave amplitude a'. In case of small wave amplitudes, 
according to the linear theory, we have the "natural" angular frequency a^ ~ a/*? |k'|. 
Then, the above resonance condition becomes 



a 



(J, 



0' 



I.e. 



E 

n=l 



m„. ar, 



cr. 



0) 



J2ml>l. 



(90) 



(91) 



n=l 
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Physically speaking, the wave resonance occurs when the nonlinear-interaction angu- 



lar frequency a' 



Y. ™n (^r. 



n=l 



of the corresponding nonlinear-interaction wave with 



wavenumber k' = ^ m„ k„ equals to its "natural" angular frequency ctq. Note that, 



n=l 



different from the nonlinear-interaction angular frequency a' that is a kind of sum of 
angular frequencies of primary waves, the "natural" angular frequency ctq of the corre- 
sponding wave number k' depends only upon the wavenumber k' and its amplitude a', 
but has nothing to do with the angular frequencies of primary waves. Thus, in general, 
the nonlinear-interaction angular frequency a' is not equal to the "natural" angular 
frequency o"q of the nonlinear-interaction wave with wavenumber k'. So, the wave 
resonance condition is indeed rather special. This physical explanation agrees well 
with the traditional resonance theory. So, fl^U]) and ( MT\i reveal the physical essence 
of the gravity wave resonance. 

Although (191 p is derived from the resonance condition (ISO]) for small wave ampli- 
tudes, this physical mechanism of gravity wave resonance has general meanings and 
holds for large wave amplitudes even if cr„ ~ \fgKi is not a good approximation. So, 
f l9T]) is also the wave resonance condition for arbitrary number of primary waves with 
large amplitudes. It should be emphasized that the wave resonance condition f l?!]) log- 
ically contains the resonance condition f lSU]) for arbitrary number of small-amplitude 
waves and Phillips' resonance condition ([1]) for four small-amplitude waves. Thus, it 
is rather general. 

When the wave resonance condition f l9T|) is satisfied and the wave energy transfers 
from the primary waves to a resonant one, the amplitudes of primary waves decreases 
and the amplitude of the resonant wave increases. Therefore, the angular frequencies 
(j„ of each primary waves decrease but the "natural" angular frequency of the res- 
onant wave increases so that the resonance condition (lOTjl does not hold any more. 
As a result, the "natural" frequency ctq departs more and more from the nonlinear- 
interaction frequency a', and the nonlinear-interaction wave gains less and less energy 
from the primary waves, until the whole wave system is in equilibrium. This explains 
why a resonant gravity wave has finite value of amplitude. As mentioned above, a 
resonant simple pendulum acted by an external force with the phase difference vr, 
as shown in FigJTT| loses its energy so that the maximum angel of oscillation 9max 
decreases. Similarly, when the wave resonance condition f l9T]) is satisfied, it is also 
possible that the wave energy transfers from the resonant wave to primary ones so 
that the amplitude of resonant wave decreases and the amplitudes of primary waves 
increase: this explains why the amplitude of a resonant wave may be much smaller 
than those of primary ones, as shown in Table | 



5 Concluding remark and discussions 

The main findings and concluding remarks are outlined below. 
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First of all, based on a analytic technique for strongly nonlinear problems, namely 
the homotopy analysis method (HAM), a multiple- variable technique is proposed and 
applied to give convergent series solution of a fully developed system of arbitrary 
number of primary periodic traveling waves. Different from perturbation techniques 
used by Phillips [Ij and Longuet-Higgins [2] , this multiple- variable technique does not 
depend upon any small physical parameters, and besides provides a convenient way 
to ensure the fast convergence of solution series. By means of this multiple-variable 
technique, the time t does not explicitly appear for a fully developed wave system: this 
not only greatly simplifies solving the problem mathematically, but also contributes a 
lot to revealing the physical meanings clearly (some users of the HAM solved nonlinear 
wave-type PDEs by simply expanding the solution in Taylor series with respect to the 
time t. Unfortunately, this often leads to very complicated solution expressions with 
rather little physical meanings). Thus, the homotopy multiple- variable method has 
general meanings and can be widely applied to different types of nonlinear problems 
in science and engineering. For example, although this method is used here for fully 
developed gravity waves, it can be applied to study the evolution of nonlinear waves far 
from equilibrium (such as the famous natural phenomena about freak wave |19H2T]), 
as long as we introduce a "slow" time-scale t = t/{l + t) as an additional variable to 
describe the time-dependent variation of the wave amplitude and angular frequency. 

Secondly, by means of the homotopy multiple- variable method, we illustrated that 
the amplitudes of all wave components of a fully developed wave system are finite 
constants, even if the resonance condition is exactly satisfied. Besides, we revealed, 
maybe for the first time, that a fully developed resonant wave system may have 
multiple solutions. Especially, it is found that the amplitude of resonant waves might 
be much smaller than that of primary waves, and that a resonant wave may contain 
only the few of the whole wave energy. These results differ from some of our traditional 
thoughts, but strongly suggest that, due to nonlinear interaction, the evolution of a 
multiple-wave system might be rather complicated. At the end of §4, some physical 
explanations for these results are given. 

Third, by means of the homotopy multiple-variable method, we derived two gen- 
eral wave resonance conditions fl5U|) and (jHI]) for arbitrary number of primary periodic 
traveling waves: the former holds for small-amplitude waves, but the latter works even 
for large-amplitude waves. These two resonance conditions logically contain Phillips' 
resonance condition ([1]) for four small-amplitude waves, and thus are more general. 
Especially, the wave resonance condition ( 19T|) opens a new way to study the strongly 
nonlinear interactions of more than four primary traveling waves with large ampli- 
tudes. 

Mathematically, our computations suggest that, for a fully developed wave system 
with small amplitudes, the main wave energy distribute in the wave components 
whose eigenvalues to the linear operator (1761) or (187|) are zero (or close to zero). 
Physically speaking, the primary and resonant waves contain the main of the wave 
energy. However, as the wave amplitudes increases so that the nonlinearity becomes 
stronger, the primary and resonant waves as a whole contain less and less percentage 
of wave energy, as shown in Fig. HOl 



39 



There are some open questions. The resonance condition (19T|) for arbitrary num- 
ber of travehng waves with large amphtudes is given from the physical view-points of 
resonance. Although this resonance condition explains very well why the amplitude 
of a resonant wave is finite and why it can be much smaller than those of primary 
ones, it should be verified by experiments or other analytical/numerical approaches. 
Besides, it is worthwhile studying the evolution of a system of multiple traveling waves 
far from equilibrium. 

Finally, it should be pointed out once again that the homotopy multiple- variable 
method proposed in this article is more general than the famous multiple-scales tech- 
niques in perturbation theory. By means of the perturbation multiple-scale technique, 
one often rewrites a unknown function f{t) in the form F{Tq, Ti, T2), where 

To = t, T, = et, T2 = eH 

denote different timescales with the small physical parameter e. By means of this 
traditional multiple-scale technique, a nonlinear problem is often transformed into a 
sequence of linear perturbed problems via the small physical parameter e. Using the 
homotopy multiple-variable method, we can also rewrite /(t) by F{E,qi ^11^,2) with the 
definition 

in = e" t. 

However, different from the multiple-scale perturbation techniques, we now do not 
need any small physical parameters to transform the original nonlinear problem into 
a sequence of linear sub-problems. Furthermore, it is easy to get high-order approx- 
imation by our approach, as illustrated in this article. Especially, if the multiple- 
variables are properly defined with clear physical meanings, this method is helpful to 
get results with important physical meanings. This work illustrates that the homo- 
topy multiple-variable method can overcome the restrictions of traditional analytic 
methods and besides it belongs to the times of computer. It seems that the homotopy 
multiple-variable technique can be applied widely to solve different types of strongly 
nonlinear problems in science and engineering. 
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Appendix A 
The detailed derivation of dH]) and d43b 



Write 

with the definition 
Then, 



E"-«' =E 



l-^m,n Q 



i=l 



/il,n(6,6) =^n(6,6), n>l. 



hoo \ "^+1 / +00 

E".^' =(E 

J=i / 

hoo / s — 1 



hoo 



E' 



\n=m / \j=l 

+00 /n — 1 

X] <?'( X] /^™." ^^-" I = X] "?" I X] '"'"•^ '^"-^ 

s=m+l \n=m / n=m+l \i=m 

+00 

/ ^ ^m+l,n Q , 



which gives 



n-l 
/^m+l,n(^l,6) = X]/^m,i(^l,6) ^n-i(6,6), m > 1, n > m + 1, 



42 



(92) 
(93) 



(94) 



I.e. 



n — 1 



/im,n(6,6)= ^ /im-l,i(^l,6) '7n-i(^l,6), m > 2, U > 171. (95) 



j=m— 1 



Thus, by means of fl93|) and fl95|) . one can easily get /i^.^^ even for large m and n. 
Define 



f)i+j / 1 /5m 



By Taylor series, we have for any z that 

1 d'^th 



z=0 



-1-00 ^ 
m=0 ^ 



ml dz^ 



z=Q 



z 






m m. 



m=0 



and 



a^+^( 



E 



^i+i / 1 ff 



^ei^e^' fr^o^^i^^2 v"^! 9z 



2=0 



m.=0 



m m 



(96) 



(97) 
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Then, on z = r)(^i,^2; q), we have using ( 192|) that 
gi+j 1 +°° /+°° 



where 



I.e. 



dmi 



m=0 \s=l 

+00 / + 

m=l \s=m 

+00 / s 

/ n,0 I \ ^ s I \ ^ / n,m 

s=l \r?i=l 

+00 / rn 

I n,0 I \ ^ m I \ ^ ; n,s 

"Pi J +2^^ z^V'i,; /is,. 

m=l \'5= 

+00 



,s=l 



m=0 



n,0 



/?"'.' 
r-^l.1 



A 



,n,0 



m 



"^i^lf ^^s,m, m > 1, 



s=l 



n,l 
n,2 
n,3 



I n,l I r> / n,2 I 3 / n,3 



A 
A 
A 

on, 4 ; n,l I /o I 2\ ; n,2 , o 2 ; n,3 , 4 ; 

A,; = ^4 V^i,; + (2r7ir73 + r/2 j V^^j- + 3r/ir/2 V^^j- + r/^ V^^ 



n, 



4 „;,n,4 



(98) 

(99) 
(100) 



Similarly, on z = fi{C,i,C,2', q), it holds 



gi+j 


\dz J 


y^ ^i+i / 1 5'»+l0„ 




deM 


;^o^a5e^ U' 5^"^+' 


2 = 0. 






+00 /+00 
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m=0 
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and 



where 



I.e. 



8'+^ {d'^cj). 



diidii \ dz 



7. 



n,l 
ra,2 

7j,i 

71, 4 



7. 



^"'0 
'i,j 



7 



«j 



2^1 Cf , 



«j 



+2^ 5*+J / 1 9™+20r 

+ 00 
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^(m + l)(m + 2) V" 

m=0 

+ 00 
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^: 
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m 
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and 
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(105) 
(106) 



Note that the explicit expressions of /3j"}™(G) G), 7j"j™(^i, G), '^"^'"(G) G) defined above 
can be easily obtained by symbohc software such as Mathematica, Maple and so on. 
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Then, on z = ?)(Ci,^2; <?), it holds using (!98|) that 



71=0 

+00 +00 

n=0 m=0 
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m+n 
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where 
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Similarly, we have 
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Similarly, on z = 57(^1, ^2; ?), we have 
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Then, on z = r)(^i,^2; <?), it holds using fllOQp and (II 111) that 
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where 
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Besides, on z = fj{^i,^2] q), we have by means of (11091) . f illip and (11121) that 
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Furthermore, using (IT09|) . dTTT]) . (ITT9|) and (IT2TD . we have 
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Then, using ( HW^ . f lTTT]) . dTTTj) . f lTTg]) and f lT^ . we have on z = 77(6,6; ?) that 
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respectively. Then, on z = ^(^1)^2;?), it holds due to the linear property of the 

operator (!28|) that 

+00 

>C(0-0o)=$^5n(ei,6)g", (129) 

n=l 
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where 

_ J 0, whenn < 1, 
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Substituting f ll3ip . f ll25p into ( 129|) and equating the like-power of q, we have the 
boundary condition: 

Srai^i, 6) - Xin 5'„_i(6, 6) = cq A* „i(^i, ^2), m > 1. (133) 
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5^1 56^6 5^2 ^^ 

Substituting the above expression into ( 1133^ gives the boundary condition on z = 0: 

>C (0m) = Co A^_i + Xm 5'm-i - 5"™, m > 1, (136) 

where £ is defined by ( H5l) . 

Substituting the series (1371) . (I109p and (IllSp into (13T1) . equating the like-power of 
g, we have 

'7m(^i, 6) = Co A^„i + Xm Vm-u m > 1, (137) 

where 

Al = Vm-^ (^1 0^° + cr2 0^^ - rm,o) . 
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Appendix B 
A brief proof of the Convergence Theorem 



Proof: The potential function is a linear combination of the eigenf unctions 

\l/m,n = exp(|mki + nk2\z) sin(m,^i + 12^2) 

for integers m and n. Thus, (p automatically satisfies the linear Laplace equation flT6|) 
and the boundary condition fl^Tl) . 
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Setting g = 1 in flTOTl) and fiTUgi) gives on z = r/(^i, ^2; 1) that 



E^i 







96 
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m=0 m=0 

Similarly, setting g = 1 in flllSp gives on 2; = 17(^1,^2; 1) that 

+00 



db 
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m,o = 2 V0(6, 6, z; 1) ■ V0(ei, 6, ;2; 1). 



Thus, we have on z = r)(^i, ^2; 1) that 
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^1 77Z + fT2 — ttW^i, 6, z; 1 ■ V0^i, G, z; 1 
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Besides, when ^ A^ = 0, setting g = 1 in (11251) gives 

m=0 



+00 



^At = Ar[0(6,6,^;i)]=o 



(139) 



m=0 



on the free surface z = 17(^1, ^2; 1)- Note that, according to fl36l) and fl37|) . 0(^i, ^2, z; 1) 
and 57(^1,^2;!) denote the homotopy-series fl38l) and fl39|) . respectively. Note also 
that (I138p and (I139p are exactly the two boundary conditions on the free surface. 
Therefore, the homotopy-series fl55]) and fl5I?l) satisfy the original governing equation 
(TT6|) and all boundary conditions (fT9|) . (120|) and (12T!) . This ends the proof. 



